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We present fully nonlinear dissipative fluid dynamics simulations of a trapped two-dimensional 
Fermi gas at unitarity using a Lattice Boltzmann algorithm. We are able to simulate non-harmonic 
trapping potentials, temperature-dependent viscosities as well as a discretized version of the ballistic 
(non-interacting) behavior. Our approach lends itself to direct comparison with experimental data, 
opening up the possibility of a precision determination of transport coefficients in the unitary Fermi 
gas. Furthermore, we predict the presence of a non-hydrodynamic component in the quadrupole 
mode, which should be observable experimentally. 
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I. INTRODUCTION 


Cold Fermi gases at unitarity are examples of so-called ’Strongly Interacting Quantum Fluids’ (SIQFs), which 
share the common characteristic that they flow around obstacles almost without friction. Other experimental ex¬ 
amples of SIQFs seemingly include hi^ temperature superconductors [ij and hot quark gluon plasmas generated in 
ultrarelativistic heavy-ion collisions [JS] , while the cleanest theoretical example seems to be given by black holes . 

The last decade has seen an unprecedented development in precision experiments in these ultracold quantum gases 
at unitarity, both in two and three dimensions, i.e. |5l-lll|. So far, theoretical descriptions of the collective behavior 
observed in these experiments has mostly been based on (analytic) linear or scaling solutions to hydrodynamic or 
kinetic equations, see e.g. Refs. [l^ - [l^ . While these solutions have the advantage of being fully analytical, they 
also have the drawback of applying only to specific, idealized situations, such as harmonic trapping potentials, fixed 
temperature dependence of transport coefficients and hydrodynamic behavior in the cloud’s low density corona. 

In the present work we follow a completely different approach in that we intend to set up numerical ‘experiments’ 
that are capable of simulating the behavior of cold atomic gases for setups that include, but are not limited to, the 
experimentally realized scenarios. We expect that by being able to simulate realistic (e.g. non-idealized) configurations 
as in experiment our method could lead to much more precise extraction of transport properties from experimental data 
than using analytic solutions. Moreover, by being able to simulate configurations that are not (yet) experimentally 
realized we can hope to point out and possibly predict interesting cloud behavior. This work is meant to be a first 
step in this direction, in that we limit ourselves to discuss only the case of the collective modes of a trapped, two- 
dimensional Fermi gas at unitarity and ignoring heat conduction effects. Generalizations of our approach to three 
dimensions, Bose gases and inclusion of heat conduction are straightforward and will be considered in future work. 

The simulations we perform are meant to describe the bulk evolution and transport phenomena occurring in an 
atomic cloud or clouds, in particular ignoring the effect of quantum transitions and quantum tunneling between 
spatially separated clouds, and outside the superfluid regime. In the absence of dissipation, such dynamics would be 
well described by one-component hydrodynamics. Given that we are interested in the effects of dissipation (transport), 
and that dissipation is always important in the low density corona of the cloud or clouds, a description in terms of 
the so-called Lattice Boltzmann (LB) algorithm [l^ may be well suited to offer a realistic description of the atomic 
cloud dynamics. Despite their name and origin, it should be pointed out that Lattice Boltzmann simulations are not 
limited to situations with well-defined quasi-particles or weak coupling situations, which is why they are well suited 
for simulating SIQFs. 

Compared to recent results in the literature, our study is closely related to Ref. [l^ , where the Boltzmann equation 
was solved numerically using the test particle method. Compared to Ref. [l^ , our method has the disadvantage of not 
accurately describing the non-interacting regime of the cold atom gas on a quantitative level. However, the present 
approach has the advantage of allowing for arbitrary equations of state, the straightforward simulation of shock waves 
as well as being considerably cheaper in terms of computational cost. 

Furthermore, the present work is related to Ref. [13, where non-hydrodynamic elements from the Boltzmann 
equation were used to improve the hydrodynamic evolution equations in the non-interacting regime. Similar to 
Ref. [T^, the algorithm used in Ref. [T3 is superior to our approach in that the non-interacting regime is treated 
exactly, but unlike our setup does not allow for non-ideal equations of state. 

Moreover, our work is similar to Refs. MM where the Boltzmann equation was solved numerically in the relaxation 
time approximation. Compared to our work, the method in Refs. [iSl Il9l | can obtain reliable results also in the non¬ 
interacting regime, but our approach is computationally cheaper. 

Finally, our work is closely related to Refs. [13,113 where the Boltzmann equation was studied in the moment 
approximation. This is because the lattice Boltzmann method is essentially a moment approximation to the Boltzmann 
equation packaged in a computationally highly efficient form. However, compared to the moment approximation of 
the Boltzmann equation, the lattice Boltzmann framework offers more flexibility by allowing for non-ideal equations 
of state as well as interaction terms that are not realizable in a particle picture. In this sense, the lattice Boltzmann 
framework is closer to hydrodynamics than to the actual Boltzmann equation, being an effective theory for low 
frequency, small wave-number transport. 

This work is organized as follows: In section|H]we give a derivation of the Lattice Boltzmann framework for the use 
of simulating cold trapped quantum gases at unitarity. Section IHII contains our results for trapped two-dimensional 
gases with ideal and non-ideal equations of state in harmonic and Gaussian traps. We present our conclusions in 
section m and include details about our numerical scheme in an appendix. 
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II. THE LATTICE BOLTZMANN FRAMEWORK 


A. Kinetic Description of Trapped Cloud 


The Boltzmann equation for a single particle distribution function /(t,x, v) is given by 


a* + V • V 


—yu(x) ■ 

m 


f = c[f] 


( 1 ) 


where v,m are the particle’s velocity and mass, U{x) is the trapping potential and C[f] is the collision term that 
depends on the particle interactions. From the single particle distribution function one can define the local mass 
density p, macroscopic (fluid) velocity u and total (kinetic plus internal) energy density e through the following 
integral moments: 


p = m 


d^vf , pu = m 




m 



d^vv^f, 


( 2 ) 


where D is the number of space dimensions. 

It is well known that if the interactions are strong enough, then the dynamical evolution of the macroscopic quantities 
p, u, e will obey the equations of (Navier-Stokes) hydrodynamics. However, in this limit the microscopic details of the 
particle interactions become unimportant, and the hydrodynamic evolution for p, u, e stay unchanged if one replaces 
the (complicated) collision term with a simple BGK-type ansatz as long as this ansatz obeys the conservation of mass, 
momentum and energy: 

C[/] = -^^^, (3) 

tr 

where /eq is the local equilibrium distribution function and tr is the (local) relaxation time. It is straightforward to 
find static equilibrium solutions to Eq. o, which in units where c = kb = h = 1 can be expressed in terms of the 
macroscopic variables as 


, , pe'^^ f 1,D = 2 

/eqH, X, vj - = 3 ’ 


( 4 ) 


with c^(T) = ^ the local speed of sound squared and a proportionality constant that depends on the number of space 
dimensions D. Note that this implies e = ^pv? + ^pcl{T). The careful reader will at this point worry about the 
choice of a Boltzmann-type distribution function (j3]), which may seem a bad approximation for describing quantum 
gases where actual particle distribution functions should be given by Fermi-Dirac or Bose-Einstein statistics. There 
are two important points to consider for this issue. One, our description is aimed at the evolution of the macroscopic 
system variables and we do not aim for a correct description of the actual particle distribution. Two, for the bulk 
evolution of the system the information about the actual particle distribution enters only through the equation of 
state (e.g. the relation between pressure and density), and we will demonstrate how to correctly implement this 
information in the following section. 

As a consequence of the unimportance of the microscopic collision term for the hydrodynamic evolution of the 
system, the original restriction of Eq. m to well-separated particle degrees of freedom can be lifted, since Eqns. o,© 
no longer explicitly refer to particles anymore. Thus, Eqns. ©,© can be viewed as a system of effective equations 
describing the dynamic evolution of the macroscopic system variables p, u, e in the limit of strong interactions {tr —>■ 0). 
In particular, this implies that Eqns. o,® can give an accurate effective system description even in situations where 
the original Boltzmann equation ® may no longer be well-defined. Conversely, in the case of weak interactions 
(in particular in the ballistic regime tr —>■ oo), the ansatz ©, while qualitatively correct, may offer only a poor 
quantitative approximation to the exact microscopic collision term C[f]. These considerations set the regime of 
applicability for the Lattice Boltzmann framework outlined below. 

Let us now rescale coordinates so that we work in dimensionless units adapted to the cold atoms system characterized 
by some transverse size R± and some frequency lox : 


so that m becomes 


= xR± , t = t/uj ±_, V = ■vR±uj± , u = , 
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In the following, the initial condition considered is that for a static cloud in equilibrium, which to a good approxi¬ 
mation will be isothermal with an initial temperature Tq. Introducing in addition a calculational (constant) parameter 
cli which will be referred to as the ’lattice speed’ below, we chose our unit system parameter R± to be given as 


R± 



( 6 ) 


such that 


di + v-V 


^VU{x) ■ V 
Jo 


(v) 


/-/eg 

TrUJ± 


(7) 


with 9 = ^. Note that in these units, for a cloud that is at rest with temperature Tq and trapping potential U, the 
mass density initially is given by 


p oc exp 



( 8 ) 


B. From Boltzmann to Lattice Boltzmann 

So far we have discussed a treatment of trapped atomic gases in continuum kinetic theory. However, it turns out 
that as long as we are only concerned with the evolution of the macroscopic variables ([2]) we do not actually need to 
keep the full continuum information of /(t,x, v) in velocity space. Specifically, since ([2]) refer only to a finite number 
of moments of the distribution function, we can obtain an exact representation of these integrals by expanding / 
in a series of polynomials which are orthogonal and complete with respect to some reference distribution function 

.y.2 

(see e.g. [l^ for details of the Lattice Boltzmann method). In particular, choosing e as a reference distribution 
function, the orthogonal polynomials are just the well-known Hermite polynomials. In compact index-notation, we 
collect these orthogonal polynomials into tensors - where n is the order of the polynomial and the tensorial 
indices ii, i 2 , ■ ■ ■ run from I to the number of space dimensions. For instance, we have 

Po{v) = 1, Pl{v) = P , P 2 ^{v) = (9) 

Because the basis functions are polynomials, one can use quadrature-type rules to exactly represent the integrals ([2]). 
The nodes of polynomials select specific, optimized values for the values of the velocity vectors v, such that one is 
left with a velocity lattice rather than a continuous collection of velocities. Labeling the individual velocity vectors of 
the lattice through the label s (running from 1 to Q, the total number of velocity vectors in the lattice), the relevant 
moments become 




pu 
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( 10 ) 




where Wg are suitably chosen integration weights and we have introduced the pressure P{n,T). It is worth stressing 
that while the continuum solution /eq in (l7|) demands an ideal gas equation of state P{n, T) = nT, this restriction is 
(to some degree) lifted in the discretized version /s(t,x). In particular, this implies that non-ideal equations of state 
P{n,T) can be simulated using the LB method. 

In the following we will work with a previously defined velocity lattice that is known as D2Q25 (D=2 space 
dimensions, and Q = 25 velocity vectors), for which = 1 — The individual velocities and weights for this 
lattice are given in Table IIIBI 
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TABLE I. Left: Velocities and weights of the D2Q25 lattice, where to = ^(4 + VTO), ti = ^(8 — V^di), ts = — SVTO). 

Right: Illustration of the velocity vectors for the D2Q25 lattice. Note that for the D2Q25 lattice, — 1 — 


Expanding the equilibrium distribution function /eq in terms of the reference distribution function and up to third 
order in polynomials (cf. Ref. [1^) gives 
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( 11 ) 


with a proportionality constant that is specified in Eq. Q. It is straightforward to verify that truncation of /eq at 
this order of expansion still leads to exact hydrodynamic evolution equations for the macroscopic quantities p, u, e 
(for small gradients of the temperature). Eor a trapped gas, we also need an expansion of the force term in Eq.([7]) in 
terms of orthogonal polynomials: 


-2 

F • = e~^Y. x), 
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F = |iVC/(x), 


( 12 ) 


where we have suppressed the tensorial indices for simplicity. Using the orthogonality of the polynomials, a straight¬ 
forward calculation gives the coefficients a„ and finally the representation of the force term as 


F • = - 
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(13) 


With the equilibrium distribution and the force term suitably discretized, one now needs to specify the discretization 
of the space and time derivatives in Eq. d?]) . Here we use the simplest version which is to rewrite 


[a^ + v-v] /(t,x)~ 


f{t + St,x + vSt) - f{t,x) 
St 


(14) 


Since the momentum lattice D2Q25 is space-filling, this particular choice implies that if we discretize space on a cubic 
lattice, at every time step ’particles’ stream from one lattice site to the next lattice site, so there is no need for any 
interpolation schemes for Eq. m- However, note that (1141) is only exact up to first order in derivatives. Taking into 
account second-order derivatives finally leads to an evolution equation of the form 


fit + 51 X + v5t) = fit, x) (l - h) + /eq(/, x)H + 5t¥ ■ , 


(15) 


with a force F, relaxation term H as well as macroscopic variables that include numerical modifications to continuum 
expressions. Denoting these (numerically corrected) quantities by a tilde, one finds 


n = n, 



= 


1 

2 


St 



T = T. 


(16) 


Note that for the evolution of Eq. dTSl) . one first calculates e.g. nu = J2s'^^'Wsfs, and then performs the numerical 
correction, before using u to construct feq through Eq. (El- 

Given a force F and a temperature/density dependent relaxation time the simulation algorithm then can be 

summarized as follows: 
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1. Select an initial condition for / (e.g. from Eo. lfTTIl i 

2. Free-stream the components of / to neighboring lattice sites through f{t + St,x + vSt) = fit, x) 

3. Calculate macroscopic variables n, u, T for this new configuration 

4. Calculate /eq, F from the macroscopic variables 

5. Correct f{i+6t) according to Ea. dTK]) 

6. Repeat from Step 2 

For completeness, we also note that the relaxation time tji is controlling the simulated ratio (and tempera¬ 
ture/density dependence) of the shear viscosity rj over the pressure P of the cloud through the relation (cf. Ref. M) 

rR = r,lP- ( 17 ) 

Note that other transport processes (e.g. spin diffusion) would have transport times that differ from and also the 
relation between diffusion constant and diffusion transport time would not be identical to Eq. ini- Here we limit 
ourselves to just momentum transport controlled by the shear viscosity. 


III. RESULTS 


A. 2D Ideal Gas in an Harmonic Trap: Analytics 

Solving Eq. © in two dimensions for a static equilibrium solution for / in an arbitrary trapping potential with 
T = Tq and an ideal gas equation of state P{n,T) = nT gives 


. U(x) 


/o(x,v) = 


]{T^) 


(18) 


The simplest case to study is that of a symmetric harmonic trapping potential U{x) = A simple solution to 

Eq. (O with harmonic potential is the case of center-of-mass oscillations (sloshing mode), which can be written as 


/(t, X, v) = /o {xi - c(t), Vi - c'{t)) , c{t) + c"{t) = 0 . 

A general ansatz for the sloshing mode equation of motion is 

c{t) = cos{iost + const), 

and in the idealized case at hand one finds Eg = 0,a;s = 1 (in units of the base frequency w_l). 
A different solution to Eq. 0 is given by a scaling ansatz [2^ , 


/(t,x,v) = /o 


Xi Vi-Xib'^{t)/h{t) 


hit)' 






0'(t) + 2^h(i) = -^’ ^ 


h(i) 


rRU!± 


(19) 


( 20 ) 


( 21 ) 


where 9 = -^ Si=i (note that this solution is valid also for D ^ 2). In the limit of small perturbations 6^(7) = 
1 + Shit), one can separate the equations into a breathing mode SBit) = ^ quadrupole mode 

SQit) = ^ Usijqg furthermore initial conditions with i56((0) = 0 one finds 


SB"it) + ASBit) = 0, SQ"it) + 2SQit) + truj± iSQ'"it) + ASQ'if)) = 0 . 

From these linear equations, one clearly can identify an undamped (Fs =0) breathing mode oscillation 

SBit) = /3e“^®‘cos(iCBt + const) 


( 22 ) 


(23) 


with frequency = 2 (in units of the base frequency u;j_) and a quadrupole mode with a frequency wq that is 
in the hydrodynamic limit tru;± —^ 0, while it increases to the same frequency as the breathing mode in the free 
streaming limit tru}± —>■ oo. A fully analytic solution to Eg. (1221) for constant uj±tr is straightforward, and given by 


SQit) = ye cosiwQt + const) + Se 


(24) 
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FIG. 1. Frequency and damping rates for the sloshing, breathing and quadrupole mode c,SB,5Qj respectively, of an ideal gas 
of atoms in an harmonic trap. Shown are exact analytic results in the limit of small amplitudes (lines) and fitted values from 
a fully numerical LB simulation (symbols with systematic errors from the infinite volume extrapolation). The ’hydrodynamic’ 
regime is the limit of oj±tii —>■ 0 (e.g. linear rise of Fq^q with ljJ±tr), whereas the ’ballistic’ (or ’free streaming’) regime 
corresponds to ljJ±tr —^ 00 . The plot suggests a turnover point from the hydrodynamic regime to the ballistic regime at around 
0 J±rR ~ 0.5. The comparison between analytic and numerical LB results indicates that the fully numerical simulation reproduces 
the analytic results rather well even for values of the relaxation time that are approaching the ballistic (free streaming) limit. 
The only exceptions to this agreement are the extraction of the non-hydrodynamic quadrupole damping rate r( 5 p, as well as 
the quadrupole frequency wq for ljJ±tr > 0.5, which are qualitatively similar, but quantitatively different in the analytic and 
numerical results (see text for details). 


where the explicit form for Fq 0 , Fq 1 , wq is lengthy. A plot of these quantities as a function of is given in FigU] 

and we note that the asymptotic behavior is 


< 1 : Fq_o ~ uj±tr , 
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(25) 


Inspecting Eq. (I24|) . it becomes clear that it is a superposition of two different modes: a well-known ’hydrodynamic 
quadrupole mode’ which is becoming dominant in the hydrodynamic limit uj±tr —>■ 0 and a ’non-hydrodynamic’ purely 
damped mode. The non-hydrodynamic mode is a feature that is common to evolution equations beyond Navier-Stokes 
(such as second-order hydrodynamics, cf. 1^). 

Because the derivation of Eq.(l7|) was carried out close to the hydrodynamic limit, we do not expect the analytically 
calculated value of FQ_i(w_Lr/{) to be quantitatively reflected in any experimental measurement. Similarly, because 
the LB framework truncates the continuum Boltzmann equation onto a finite number of basis functions, we do not 
expect the numerical LB result to match the analytic value for Fn 1 , either. However, because the derivation is still 
qualitatively sound in this regime, we expect a term such as to also be present and observable in both the 

numerical LB simulation and experiments of trapped atomic clouds. Furthermore, we argue that extraction of the 
coefficient FQp(a;_LT/j) from data could be very interesting because it is this coefficient which will indicate the radius 
of convergence of the hydrodynamic approximation. This is evident from recent progress in the context of relativistic 
fluid dynamics, where a term such as with non-hydrodynamic dependence on the relaxation time (1251) 

indicates the presence of so-called quasi-normal mode behavior [261 l27l| (while there is only one such term in Eq. (1241) , 
in practice we expect experimental signals to contain an infinite series of terms of the form possibly also with 

oscillating components, with Fq^„ oc n for n ;g> 1 originating from an infinite tower of quasinormal modes). Thus, by 
experimentally determining Fq i (or possibly also information about higher order quasinormal modes) one can expect 
to learn about non-hydrodynamic behavior in strongly-coupled quantum fluids and we encourage experimentalists to 
consider this option in future work. 


B. 2D Ideal Gas in an Harmonic Trap: Nnmerical Simulation 

In the previous subsection we have considered the analytic solution to the breathing and quadrupole modes of a 
two-dimensional Fermi gas in an harmonic trap. Here we proceed to simulate this setup using the Lattice Boltzmann 
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FIG. 2. Two-dimensional snapshots of the cloud density profile n(x, y) (normalized by the central density no) for the breathing 
and quadrupole mode simulations. Time stamps are calculated assuming a transverse trapping frequency of (jJ± = 27r x 125 Hz 

0 - 


algorithm outlined in section lll Bl As advertised, we discretize the two-dimensional space on a square grid with lattice 
spacing St such that x = i St with i = 1... N where N is the number of gridpoints we simulate along one dimension. 
Thus, the infinite volume limit and continuum limit of our simulation correspond to taking iV —>■ oo with Si = const 
and Si ^ 0 with NSi = const, respectively. In the results we present in the following, we have performed multiple 
simulations for different volumes and resolutions; while the infinite volume and continuum limit can never be reached 
in practice, we have striven to obtain robust extrapolation of our results to these limits, and we report the residual 
error from the extrapolation procedure in our LB results (see appendix lAl for details about the numerical procedure). 

We first initialize the simulation in a configuration corresponding to the analytically tractable case presented in 
section fill A1 that is, a harmonic trapping potential with small oscillations around the equilibrium configuration. To 
this end, we define the amplitude bi(t) as the distance from the center where the density has dropped by a factor 
e“^, and we initialize the particle distribution function according to (I18I21I) with Sbi{0) ~ 10“^ or smaller. We then 
track the time evolution of the combinations SB{i),SQ{i) and perform a fit of this evolution using the analytically 
derived model functions (I23I24I) . obtaining best-fit values for the coefficients wb, wq, Tb, Lq^q, as a function of 
w_ltb in the process. The extracted values for these coefficients are displayed in Fig[T] along with the analytically 
calculated results. As can be seen from this comparison, the numerical simulation is in good quantitative agreement 
with the analytic results for the coefficients uib, Lb, Fq^Oi even for values of w_ltb which are outside the hydrodynamic 
regime. The quadrupole mode frequency wq agrees very well with the analytic result in the hydrodynamic limit, but 
starts to differ noticeably for values of u}±Tii > 0.5. We attribute this disagreement to the discretization procedure of 
continuous velocities onto the D2Q25 grid, which implies that the solution to the LB equation only corresponds to the 
solution of the continuum Boltzmann equation in the hydrodynamic limit. We plan to test this hypothesis in future 
work by employing discretization grids with a larger number of velocity vectors, which should lead to wq results that 
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are closer to the analytic values. 

As anticipated, the (non-hydrodynamic mode) damping rate Fq i is found to be qualitatively similar to the analytic 
result, but in clear quantitative disagreement. Also, note that since Fq^i becomes very large in the hydrodynamic 
regime compared to the hydrodynamic damping Fq o, it becomes more and more difficult to extract Fg i from the 
numerical simulation in the limit of a;_Lr/j ^ 1. We expect this to happen also in an experimental setup. Ffowever, 
by inspecting the analytic solutions given in Eqns. (|23)) . (l24ll . one could try to engineer initial conditions for the cloud 
that would maximize the amplitude 6, thus presumably leading to a better signal to noise ratio for extracting Fq i in 
the hydrodynamic limit. We intend to pursue this direction in a follow-up study. 

To summarize, our numerical LB simulation is able to accurately reproduce the collective behavior of a cold atomic 
gas cloud in the limit of small amplitude oscillations. This should be considered a successful test of the method. In 
the following, we will now use the numerical LB simulation to study the bulk evolution of a cold atomic gas cloud for 
situations where an analytic treatment is either not possible or difficult. 


Interacting 2D Fermi Gas in a Gaussian Trap: Numerical Simulation 


In actual experimental setups, the trapping potential is usually not harmonic. In the case of 2D Fermi Gases, it is 
more accurately described by a Gaussian potential (cf. @) 

t/(x) = I/o , (26) 

where the potential depth is Vq and the parameter a is related to the laser beam waist. If this potential is meant to 
approximate an harmonic trap close to the center x ~ 0, then cr^ = Fixing a in this way, the Gaussian trap 

corresponds to a one-parameter generalization of the harmonic trap, with ^ controlling the degree of anharmonicity 
(the case of the purely harmonic trap is recovered in the limit ^ ^ oo). 

For an ideal gas, the equation of state takes the form P{n, T) = riT, and the results in sections IlII Al IllT^ have been 
obtained by using this (idealized) equation of state. For an ideal (meaning non-interacting) Fermi gas, the equation 
of state is different from that of an ideal gas because the Fermi statistics imply a non-linear relation between the 
pressure and the density. Furthermore, in setups relevant for cold atom emeriments at unitarity, the equation of 
state is known to be different from both the ideal gas and ideal Fermi gas Q, particularly for the two dimensional 
case D = 2 (^ . In order to have a realistic description of the dynamics, we thus implement the interacting equation 
of state from Ref. [2^ in our simulations. The equation of state is constructed out of tabulated data for the density 
as a function of chemical potential /r and temperature T, n = n{fj,/T). Results are available for various values of 
the phy sical binding energy of the two-body bound state Eb, which is always present for an attractive 2D Fermi gas 
[291 [30| . From the density, the pressure can be calculated numerically through direct integration of the thermodynamic 

relation n = , matched to the virial expansion for small densities: 


P^nT(l + B2n), ^<- 5 . 


(27) 


With the equation of state fixed, the initial conditions for an isothermal atomic cloud are corresponding to solutions 
to the equations of hydrostatics (cf. Fq. ([5])), 


VP(fi,T) = nV/i(x) = -nVU(x). 

Thus, a solution to the hydrostatic equations, in terms of rescaled coordinates , is given by 

/r(x) = /To - t/(x), 


(28) 


(29) 


where fiQ, the chemical potential at the trap center, is sometimes referred to as the Fermi Energy. For fixed /xg and 
temperature Tq, the total particle number N is given by an integral over the number density. 


N = J (fxn{^{x,T)) . 


(30) 


In the idealized case of harmonic trapping potential and non-interacting Fermi gas, Fq. (1301) simplifies in the zero- 

2 

temperature limit as —>■ which is often used to define an idealized “Fermi Temperature” Tp = /tq of a trapped 

atomic gas as 


Tp = '\/2Nuj± . 


( 31 ) 
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FIG. 3. Frequency and damping rates for the sloshing, breathing and quadrupole mode c, SB, SQ, respectively, for an interacting 
2D Fermi gas in a Gaussian trap with ^ = 10. Shown are htted values from a fully numerical LB simulation (symbols with 
systematic errors from the infinite volume extrapolation where available; otherwise infinite volume trend is indicated by an 
arrow, see appendix El for details). For comparison, the analytic results for an ideal gas in a harmonic trap (cf. Fig. [T]). We 
find that the extracted frequencies in the interacting, Gaussian trap case show the same qualitative behavior as the idealized 
analytic result as a function of oj±tr, but are systematically lower. Furthermore, the extracted breathing mode damping rates 
FsiFs in the interacting, Gaussian trap case are non-vanishing, in contrast to the idealized analytic result. 


In experiments on a cloud of trapped atoms close to unitarity, the temperature of the system is reported in 
units of Tp as defined in Eq. m, finding ^ = 0.37 — 0.9 for an average number of atoms of TV ~ 2000 per 2D 
cloud 0. Once the temperature and the number of atoms is known, the chemical potential /ig at the center of an 
arbitrary trapping potential can be obtained by numerically inverting Eq. (1301) . which is the strategy we employ in 
the following. To be explicit, we choose to simulate a 2D cloud of TV = 2320 atoms and ^ = 0.45 in the following. 
At any instant in time, we assume the cloud to be isothermal, but we allow the temperature to fluctuate as a function 
of time T = T{t). Also, we choose an interaction strength corresponding to ^ = 1.0, which can be related to the 
experimentally reported quantity In(fcpa) as follows [^ : 


In (kpa) 


1 , \ Tq Eb 

2 [2rp To 


0.74. 


(32) 


A fully realistic implementation of the trapping potential would require precise knowledge of the laser beam waist 
parameters, gravitational effects as well as the magnetic field gradients (see e.g. the discussion in Ref. HI). While we 
aim to implement this in a follow-up study, for the present work we chose a reference value ^ = 10 that is comparable 
to experimental values in Ref. M- 

The above choices of parameters imply a central chemical potential over temperature ratio of ^ ~ 0.1. Note that 
for a non-ideal equation of state, the simulated viscosity over density ratio becomes (cf. Eq. (fT71) ) 


H , ( 33 ) 

n nip 

which is in general density and temperature dependent. Eor further reference, we note that ^ — 2 for uJpTp = 0.1 at 
/ig = 0. However, it should be emphasized that strictly speaking the viscosity of a two-dimensional fluid is ill-defined 
because of the presence of thermal fluctuations (cf. [321 l33l|). a topic which we intend to revisit in future work. 

Once the physical parameters have been specified, we simulate the collective behavior of the cold atomic cloud using 
the LB algorithm outlined in Sec lIII Bl We start with the initial condition of a slightly perturbed cloud in a Gaussian 
trap and track the time evolution of the sloshing, breathing and quadrupole modes, c(t), 5B{t), 6Q(t). 

Unlike in the case of ideal equation of state and harmonic trapping potential, we find signs of transient phenomena 
not captured by the solution structure given in Eans. (j20l23l24ll in LB simulations at finite volume and resolution. 
These transient phenomena could be due to ’overtones’ to the hydrodynamic sloshing, breathing and quadrupole 
modes, with higher frequencies and damping rates than the fundamental modes. However, our current numerical 
accuracy is not sufficient to distinguish real overtones from possible numerical artifacts at finite volume and resolution, 
so we intend to revisit this issue in a follow-up high precision study. 

In general, one finds that the presence of an anharmonic trapping potential will affect the evolution equations for 
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Frequency of Breathing and Quadrupole Modes Damping Rate of Breathing and Quadrupoie Mode 




Vo/To Vo/To 

FIG. 4. Frequency and damping rates for the breathing and quadrupole mode SB,SQ, respectively, for an interacting 2D Fermi 
gas with uj±rii = 0.1 in a Gaussian trap varying anharmonicity parameter Shown are fitted values from a fully numerical 
LB simulation (symbols with systematic errors from the infinite volume extrapolation). Dotted lines indicate the results for 
^ = 100. In the harmonic trap limit (corresponding to lim ^ —>■ oo), one essentially recovers the idealized results shown 
in Fig. [T] indicating that the non-ideal equation of state does not strongly affect either quadrupole or breathing mode in the 
LB simulations. Not shown is the non-hydrodynamic quadrupole mode damping Fq^i, because for ijJ±tu = 0.1 this mode is 
strongly damped and could not be unambiguously extracted from the simulations. 


the sloshing, breathing and quadrupole mode, resulting in a change of the frequencies tcg, wq,wb with respect to the 
harmonic trap values Eq. (l2^ . This effect has been noted before [s^. 

The extracted frequencies ws,wbtWq and damping rates Tg, T^, Tq^Q) from the LB simulation of collective 
modes of a two-dimensional non-ideal Fermi gas in a Gaussian trap are shown in Fig.[3l The behavior of the frequencies 
and damping rates are qualitatively similar to the case of the idealized, harmonic trap case, but shifted to lower values. 
Since in the experimental setup the sloshing mode frequency is used to calibrate wj,, such a frequency shift is not 
apparent in the experimental measurements 0- Moreover, one finds that in the Gaussian potential trap both the 
breathing mode and sloshing mode damping rates rB,rs are no longer consistent with zero, but found to be a small 
but non-vanishing value, similar to what has been found in experiment 0- 

Indeed, when considering the harmonic trap limit lim ^ —>■ oo, the results for all modes, sloshing, breathing and 
quadrupole mode, tend to the result for an ideal gas in an harmonic trap shown in Fig. [TJ This suggests that any 
non-ideal equation of state effects have a minor impact on the extracted frequencies and damping rates of the modes 
considered here. Thus, the frequencies and damping rates are essentially controlled by the value of uj±tb as well 
as the trap anharmonicity parameter In view of this we can attempt to compare the extracted frequencies and 
dampin g ra tes in our LB simulations to experimentally determined values, as has been done before by other authors 
(cf. [l9|, I35l437l| ). To do this, we need to match the value of uj±tb to the interaction strength. At finite interaction 
strength, it is reasonable to assume that the relaxation time tb is proportional to the inverse of the imaginary part 
of the scattering amplitude (cf. [13)0), thus 


uj±tr = K (^1 + ^ ln^(fc_Fa)^ , (34) 

with K a (density and temperature dependent) normalization factor. An average value of K can be estimated by 
matching the maximum of the quadrupole dam ping rate found in the LB simulations at uj±tb ~ 0.5 to the location 
at In(fcir'a) ~ 3 found in the experiment of Ref. [l0|, giving K ~ 0.12. The comparison between LB simulation with 
■^ = 10 of the quadrupole damping rate and frequency using Eq. pip is shown in Fig. [3 From this comparison, it 
can be seen that with the one-parameter fitting through Eq. (1341) . the overall agreement between the LB simulation 
and experimental values for the frequency and damping rate can be considered reasonable. However, note that the 
larger LB damping rate in Fig. [S] compared to the analytic result stems mainly from the strong anharmonicity effects 
encountered for ^ = 10 (cf. Fig. |4]), which may be larger than in the actual experimental setup. On the other hand, 
not included in Fig. [ 5 ] are effects of density-dependent ujtr, which are expected to increase the damping rate with 
respect to the analytic result . We intend to perform a more detailed comparison to experimental data, including 
an investigation of the above points as well as the effects of changing the temperature, in a follow-up study. 
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Damping Rates of Quadrupole Mode 
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FIG. 5. Quadrupole frequency and damping rate (normalized by the sloshing mode freqnency) measured in Ref. [a (’Exp’) 
compared to the results of the LB simulation with anharmonicity parameter ^ = 10 (’LB’) and the analytic result using the 
one-parameter fit given in Eq. (iMf. 


IV. CONCLUSIONS 

In this work, we have presented fully nonlinear dissipative fluid dynamics simulations of a trapped two-dimensional 
Fermi gas close to unitarity based on the Lattice Boltzmann algorithm. We were able to verify our simulations using 
the analytically tractable case of an ideal gas in an harmonic trap, finding excellent agreement in the fluid dynamics 
regime, as well as qualitative agreement in the ballistic (non-interacting) regime. Furthermore, we were able to 
simulate the case of non-ideal equations of state as well as non-harmonic trapping potentials, relevant to the study 
of collective modes in cold atom experiments. For convenience, we have made our simulation source code publicly 
available at [s^. 

Based on our simulations as well as analytic results, we predicted the presence of a non-hydrodynamic component 
of the quadrupole collective mode, which should be observable in experiments. 

We expect our work to be a step towards a fully realistic simulation of trapped Fermi and Bose gases in two and three 
spatial dimensions, close to (but not limited to) the unitary regime. Our study is complementary to other realistic 
simulation approaches, such as those of Refs. |16l42l| . We believe that comparing and combining these simulation 
results will open up the possibility of precision determination of transport coefficients from experiments of cold atomic 
gases, such as the shear and bulk viscosities and heat conductivity. 
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Appendix A: Details on Numerical Scheme: Infinite Volume, Continuum Limit, Conserved Quantities 

In this appendix, we give details about the precision of our numerical scheme, as well as the infinite volume and 
continuum limit. All simulation results shown here are for harmonic trapping potential and = 0.1. 

In Eig. ini we consider the time-evolution of the total particle number, total momentum and total energy in our 
simulation. All these quantities should be exactly conserved, but as is generally the case with numerical schemes, 
conservation is broken in the numerical evolution. This is not a problem as long as the violation of conserved quantities 
is small in overall magnitude and converging to zero in the infinite volume and continuum limit. Eig. |6] highlights 
the relative error for the total number density and energy conservation, and the absolute error for the momentum 
conservation for a fixed volume L = N6t = and increasingly better resolution (parametrized by increasing N 
from 251 to 2001 points). As can be seen in Eig. [BJ our scheme conserves particle number and total momentum to 
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FIG. 6. Time evolution of total particle number N{t), total momentum Px{t) {Py{t) is similar) and total energy E{t). Shown 
is relative error N{t)/N{0) — 1, E{t)/E{0) — 1 for particle number and energy and absolute error Px{t) for total momentum for 
various resolutions (parametrized by the number of gridpoints N = 251, 501,1001, 2001). 


machine precision, while energy conservation is dominated by resolution artifact. Note that in order to achieve exact 
energy conservation, one would need to include polynomials up to fourth order in Eq. m- In the case at hand, the 
energy conservation violation is not increasing in magnitude as a function of simulation time, so our simulations are 
long-time stable. Also, energy violations converge to zero quadratically with resolution 0{SP), so that our scheme is 
second-order convergent to the exact energy conservation limit. 

In Fig. [3 the continuum limit Si ^ 0 with NSi = const of the quadrupole and breathing mode oscillations are 
studied. One finds that the quadrupole mode is fairly insensitive to finite resolution artifacts, whereas the breathing 
mode converges in second-order to the continuum limit. Note that even though the total amplitude of the simulated 
breathing mode is fairly sensitive to the resolution (the reason being the numerical correction terms in Eq. (1161) ). 
neither the breathing mode frequency nor the damping rate show a strong sensitivity. This implies that reliable 
extractions of frequency and damping rate can be performed from simulations with rather coarse resolutions. 

In Fig. |S] we consider the long-term time evolution of the breathing and quadrupole mode for fixed resolution and 
various simulation volumes. As a consequence of conserving number, momentum and energy in our numerical scheme, 
our simulations remain stable essentially forever (note that for w_l = 27r x 125 Hz, ui±t = 300 corresponds to 380 
ms, much longer than typically studied in experimental setups). Interestingly, as can be seen in the rhs panel of 
Fig. [51 we find that at very late times the breathing mode amplitude starts to decay. This decay is a finite volume 
artifact, and can be removed by simulating larger volumes (see again the rhs panel of Fig. [5]), and we find that in the 
infinite volume limit {Si = const, N —> oo) the amplitude of the breathing mode is constant in an harmonic trapping 
potential, as expected from analytic results in continuum. 

For results on the frequency and damping rate shown in the main text of this article we used multiple simulations 
at different volume and different resolution. Using these results we perform extrapolations to the infinite volume 
and continuum limit. To be more specific, we first performed simulations at fixed volume NSi = const and various 
resolutions Si, and extracted damping rates and frequencies of interest for all of these simulations. To extrapolate 
to the continuum for we then performed one-parameter power-law least-square fits to the finite-resolution data. For 
example, using results for the quadrupole damping rate rQ^o(<5^ at various resolutions Si we obtain least-square fits 
of the form 


TQ^oiN,Si)=co{N)+ciiN)iSi}'^, n = 1,2,3,4,... , (Al) 

and select those values of n which have the overall smallest least-square value (“best fit”) and two more which have the 
second and third smallest least-square value (to quantify the quality of the fit and uncertainty of the extrapolation). 
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Convergence of Quadrupole and Breathing Mode Oscillations 
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FIG. 7. Continuum limit convergence of the time evolution for quadrupole mode {bx —bz) and breathing mode (bx + bz) for fixed 
volume and increasing resolution (parametrized by the number of gridpoints N along one dimension N = 251, 501,1001, 2001). 


The continuum extrapolated value for Tg o is then found by evaluating the fit function (lAll) for the best-fit value of 
n at = 0, and the uncertainty of the extrapolation is obtained by similarly evaluating the second and third-best fit 
function. The result of the procedure is shown for the case of r/jwj, =0.1 in the Ihs of Fig. O 

Once the continuum extrapolation has been performed for several volumes, we perform a min-y^ fit to the continuum 
extrapolated data including error bars. For example, we use 


TQ^o(.N-\Si = 0) = do + , 


n = 1, 2, 3,4,... , 


and select those values of n which have the overall smallest value (“best fit”) and two more which have the second 
and third smallest value (to quantify the quality of the fit and uncertainty of the extrapolation to infinite volume). 
An example of the infinite volume extrapolation of the continuum extrapolated data is shown in the rhs of Fig. 121 
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FIG. 8. Left: Long-term evolution of quadrupole mode {bx — bz) and breathing mode {bx + bz) for N = 751, St = Right: 

Infinite-volume limit for even longer time evolutions for the breathing mode for Si = and N = 601, 651, 701, 751. 
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Infinite Volume Extrapolation of Damping Rate 
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FIG. 9. Left: Continuum extrapolation of the quadrupole mode damping rate Fq^o for constant volume N6i = Right: 
infinite volume extrapolation of the continuum extrapolated damping rate values (error bars from continuum extrapolation are 
smaller than symbol sizes). 


After this extrapolation procedure we thus obtain the continuum and infinite volume extrapolated quantities such as 
lim^t_j.Q lim^v-s-oo rQ_o(A^j including uncertainty estimates from the extrapolation procedure. Wherever possible, 
we report results for these extrapolated quantities (rather than results at hnite volume or finite resolution) in the 
main text of this article. 

However, in some cases, such for anharmonic traps and non-ideal equations of state studied in the main part of 
the text, the continuum extrapolation is computationally too demanding for our present resources. In this case, we 
performed the infinite volume limit for a certain choice of parameters and then used the difference to a simulation at 
fixed volume as an indicator for the infinite volume trend at other parameter choices. For instance, Fig (TU] displays 
the volume dependence of the quadrupole damping rate for the case of a non-ideal Fermi gas in a Gaussian trap. As 
can be seen from this figure, the volume dependence of the damping rate is non-monotonic and starts to converge only 
for very large volumes. We found that the dependence of the extracted damping rate on the volume is considerably 
smaller if discarding the late-time, low amplitude simulation data. We attribute this to the fact that when the mode 
amplitude becomes low, numerical noise starts to contaminate the signal. Nevertheless, we find that for extremely 
large volumes the extracted damping rate from both procedures is equal within the statistical uncertainty. The 
difference between the extracted damping rate at low volumes {N = 500) and high volumes {N = 1700) for the case 
of TfiU! = 0.1 defines an infinite volume trend for the extracted data point at low volumes. This trend (characterized 
as both a direction and magnitude, and assumed to be independent of r^w^) has been represented as an arrow on 
the N = 500 data points for the quadrupole mode shown in Fig. |31 For completeness, we mention that the breathing 
mode does not suffer from the same severe volume dependence as found for the quadrupole mode. Thus, a standard 
infinite volume extrapolation is possible in this case. 
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Infinite Volume Extrapolation of Damping Rate 



1/N 

FIG. 10. Volume dependence of the quadrupole damping rate for a non-ideal Fermi gas in a Gaussian trap for = 0.1. 

Shown are results where the damping rate was extracted using data until late times (at which the amplitude reached one 
percent of the initial amplitude (’1% cutoff’) as well as an extraction discarding the late time date (cutting the data at an 
amplitude of ten percent (’10% cutoff’)). See text for details. 
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